function IRF=IRF_1order(MODEL,solution,StdVec,CorrMat,T_IRF)

%% Preliminaries


HX              =   solution.hx;
GX              =   solution.gx;

%% Loading Matrix on State Variables
if nargin<=3 || isempty(CorrMat)
    CorrMat         =   eye(MODEL.VarBlock.SH.Dim);
end
if nargin<=4
    T_IRF           =   41;
end
LoadMat         =   [CorrMat; ...
                     zeros(MODEL.VarBlock.XX.Dim-MODEL.VarBlock.SH.Dim, ...
                           MODEL.VarBlock.SH.Dim)] ...
                     *diag(StdVec);
%% IRF
IRF             =   struct();
for ii_SH=1:MODEL.VarBlock.SH.VarNum
    SH              =   MODEL.VarBlock.SH.VarList{ii_SH};
    IRF_XX          =   zeros(MODEL.VarBlock.XX.Dim,T_IRF);
    % IRF of XX Block
    EPS             =   zeros(MODEL.VarBlock.SH.Dim,1);
    EPS(MODEL.VarBlock.SH.LocIdx.(SH)) ...
                    =   ones(length(MODEL.VarBlock.SH.LocIdx.(SH)),1);
    IRF_XX(:,1)     =   LoadMat*EPS;
    for tt=2:T_IRF
        IRF_XX(:,tt)    =   HX*IRF_XX(:,tt-1);
    end
    IRF.(SH).XX     =   IRF_XX;
    % IRF of YY Block
    IRF_YY          =   GX*IRF_XX;
    IRF.(SH).YY     =   IRF_YY;
    
    % IRF of Individual Variables in XX and YY Blocks
    for ii_X=1:MODEL.VarBlock.XX.VarNum
        X               =   MODEL.VarBlock.XX.VarList{ii_X};
        IRF.(SH).(X)    =   IRF_XX(MODEL.VarBlock.XX.LocIdx.(X),:);
    end
    for ii_Y=1:MODEL.VarBlock.YY.VarNum
        Y               =   MODEL.VarBlock.YY.VarList{ii_Y};
        IRF.(SH).(Y)    =   IRF_YY(MODEL.VarBlock.YY.LocIdx.(Y),:);
    end
end
